home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / GELB.FOR < prev    next >
Encoding:
Text File  |  1985-11-29  |  6.9 KB  |  246 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE GELB
  5. C
  6. C        PURPOSE
  7. C           TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH A
  8. C           COEFFICIENT MATRIX OF BAND STRUCTURE.
  9. C
  10. C        USAGE
  11. C           CALL GELB(R,A,M,N,MUD,MLD,EPS,IER)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           R      - M BY N RIGHT HAND SIDE MATRIX (DESTROYED).
  15. C                    ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS.
  16. C           A      - M BY M COEFFICIENT MATRIX WITH BAND STRUCTURE
  17. C                    (DESTROYED).
  18. C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.
  19. C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.
  20. C           MUD    - THE NUMBER OF UPPER CODIAGONALS (THAT MEANS
  21. C                    CODIAGONALS ABOVE MAIN DIAGONAL).
  22. C           MLD    - THE NUMBER OF LOWER CODIAGONALS (THAT MEANS
  23. C                    CODIAGONALS BELOW MAIN DIAGONAL).
  24. C           EPS    - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
  25. C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
  26. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  27. C                    IER=0  - NO ERROR,
  28. C                    IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-
  29. C                             TERS M,MUD,MLD OR BECAUSE OF PIVOT ELEMENT
  30. C                             AT ANY ELIMINATION STEP EQUAL TO 0,
  31. C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
  32. C                             CANCE INDICATED AT ELIMINATION STEP K+1,
  33. C                             WHERE PIVOT ELEMENT WAS LESS THAN OR
  34. C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
  35. C                             ABSOLUTELY GREATEST ELEMENT OF MATRIX A.
  36. C
  37. C        REMARKS
  38. C           BAND MATRIX A IS ASSUMED TO BE STORED ROWWISE IN THE FIRST
  39. C           ME SUCCESSIVE STORAGE LOCATIONS OF TOTALLY NEEDED MA
  40. C           STORAGE LOCATIONS, WHERE
  41. C             MA=M*MC-ML*(ML+1)/2    AND    ME=MA-MU*(MU+1)/2    WITH
  42. C             MC=MIN(M,1+MUD+MLD),  ML=MC-1-MLD,  MU=MC-1-MUD.
  43. C           RIGHT HAND SIDE MATRIX R IS ASSUMED TO BE STORED COLUMNWISE
  44. C           IN N*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN SOLUTION
  45. C           MATRIX R IS STORED COLUMNWISE TOO.
  46. C           INPUT PARAMETERS M, MUD, MLD SHOULD SATISFY THE FOLLOWING
  47. C           RESTRICTIONS     MUD NOT LESS THAN ZERO
  48. C                            MLD NOT LESS THAN ZERO
  49. C                            MUD+MLD NOT GREATER THAN 2*M-2.
  50. C           NO ACTION BESIDES ERROR MESSAGE IER=-1 TAKES PLACE IF THESE
  51. C           RESTRICTIONS ARE NOT SATISFIED.
  52. C           THE PROCEDURE GIVES RESULTS IF THE RESTRICTIONS ON INPUT
  53. C           PARAMETERS ARE SATISFIED AND IF PIVOT ELEMENTS AT ALL
  54. C           ELIMINATION STEPS ARE DIFFERENT FROM 0. HOWEVER WARNING
  55. C           IER=K - IF GIVEN - INDICATES POSSIBLE LOSS OF SIGNIFICANCE.
  56. C           IN CASE OF A WELL SCALED MATRIX A AND APPROPRIATE TOLERANCE
  57. C           EPS, IER=K MAY BE INTERPRETED THAT MATRIX A HAS THE RANK K.
  58. C           NO WARNING IS GIVEN IF MATRIX A HAS NO LOWER CODIAGONAL.
  59. C
  60. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  61. C           NONE
  62. C
  63. C        METHOD
  64. C           SOLUTION IS DONE BY MEANS OF GAUSS ELIMINATION WITH
  65. C           COLUMN PIVOTING ONLY, IN ORDER TO PRESERVE BAND STRUCTURE
  66. C           IN REMAINING COEFFICIENT MATRICES.
  67. C
  68. C     ..................................................................
  69. C
  70.       SUBROUTINE GELB(R,A,M,N,MUD,MLD,EPS,IER)
  71. C
  72. C
  73.       DIMENSION R(1),A(1)
  74. C
  75. C     TEST ON WRONG INPUT PARAMETERS
  76.       IF(MLD)47,1,1
  77.     1 IF(MUD)47,2,2
  78.     2 MC=1+MLD+MUD
  79.       IF(MC+1-M-M)3,3,47
  80. C
  81. C     PREPARE INTEGER PARAMETERS
  82. C        MC=NUMBER OF COLUMNS IN MATRIX A
  83. C        MU=NUMBER OF ZEROS TO BE INSERTED IN FIRST ROW OF MATRIX A
  84. C        ML=NUMBER OF MISSING ELEMENTS IN LAST ROW OF MATRIX A
  85. C        MR=INDEX OF LAST ROW IN MATRIX A WITH MC ELEMENTS
  86. C        MZ=TOTAL NUMBER OF ZEROS TO BE INSERTED IN MATRIX A
  87. C        MA=TOTAL NUMBER OF STORAGE LOCATIONS NECESSARY FOR MATRIX A
  88. C        NM=NUMBER OF ELEMENTS IN MATRIX R
  89.     3 IF(MC-M)5,5,4
  90.     4 MC=M
  91.     5 MU=MC-MUD-1
  92.       ML=MC-MLD-1
  93.       MR=M-ML
  94.       MZ=(MU*(MU+1))/2
  95.       MA=M*MC-(ML*(ML+1))/2
  96.       NM=N*M
  97. C
  98. C     MOVE ELEMENTS BACKWARD AND SEARCH FOR ABSOLUTELY GREATEST ELEMENT
  99. C     (NOT NECESSARY IN CASE OF A MATRIX WITHOUT LOWER CODIAGONALS)
  100.       IER=0
  101.       PIV=0.
  102.       IF(MLD)14,14,6
  103.     6 JJ=MA
  104.       J=MA-MZ
  105.       KST=J
  106.       DO 9 K=1,KST
  107.       TB=A(J)
  108.       A(JJ)=TB
  109.       TB=ABS(TB)
  110.       IF(TB-PIV)8,8,7
  111.     7 PIV=TB
  112.     8 J=J-1
  113.     9 JJ=JJ-1
  114. C
  115. C     INSERT ZEROS IN FIRST MU ROWS (NOT NECESSARY IN CASE MZ=0)
  116.       IF(MZ)14,14,10
  117.    10 JJ=1
  118.       J=1+MZ
  119.       IC=1+MUD
  120.       DO 13 I=1,MU
  121.       DO 12 K=1,MC
  122.       A(JJ)=0.
  123.       IF(K-IC)11,11,12
  124.    11 A(JJ)=A(J)
  125.       J=J+1
  126.    12 JJ=JJ+1
  127.    13 IC=IC+1
  128. C
  129. C     GENERATE TEST VALUE FOR SINGULARITY
  130.    14 TOL=EPS*PIV
  131. C
  132. C
  133. C     START DECOMPOSITION LOOP
  134.       KST=1
  135.       IDST=MC
  136.       IC=MC-1
  137.       DO 38 K=1,M
  138.       IF(K-MR-1)16,16,15
  139.    15 IDST=IDST-1
  140.    16 ID=IDST
  141.       ILR=K+MLD
  142.       IF(ILR-M)18,18,17
  143.    17 ILR=M
  144.    18 II=KST
  145. C
  146. C     PIVOT SEARCH IN FIRST COLUMN (ROW INDEXES FROM I=K UP TO I=ILR)
  147.       PIV=0.
  148.       DO 22 I=K,ILR
  149.       TB=ABS(A(II))
  150.       IF(TB-PIV)20,20,19
  151.    19 PIV=TB
  152.       J=I
  153.       JJ=II
  154.    20 IF(I-MR)22,22,21
  155.    21 ID=ID-1
  156.    22 II=II+ID
  157. C
  158. C     TEST ON SINGULARITY
  159.       IF(PIV)47,47,23
  160.    23 IF(IER)26,24,26
  161.    24 IF(PIV-TOL)25,25,26
  162.    25 IER=K-1
  163.    26 PIV=1./A(JJ)
  164. C
  165. C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R
  166.       ID=J-K
  167.       DO 27 I=K,NM,M
  168.       II=I+ID
  169.       TB=PIV*R(II)
  170.       R(II)=R(I)
  171.    27 R(I)=TB
  172. C
  173. C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN COEFFICIENT MATRIX A
  174.       II=KST
  175.       J=JJ+IC
  176.       DO 28 I=JJ,J
  177.       TB=PIV*A(I)
  178.       A(I)=A(II)
  179.       A(II)=TB
  180.    28 II=II+1
  181. C
  182. C     ELEMENT REDUCTION
  183.       IF(K-ILR)29,34,34
  184.    29 ID=KST
  185.       II=K+1
  186.       MU=KST+1
  187.       MZ=KST+IC
  188.       DO 33 I=II,ILR
  189. C
  190. C     IN MATRIX A
  191.       ID=ID+MC
  192.       JJ=I-MR-1
  193.       IF(JJ)31,31,30
  194.    30 ID=ID-JJ
  195.    31 PIV=-A(ID)
  196.       J=ID+1
  197.       DO 32 JJ=MU,MZ
  198.       A(J-1)=A(J)+PIV*A(JJ)
  199.    32 J=J+1
  200.       A(J-1)=0.
  201. C
  202. C     IN MATRIX R
  203.       J=K
  204.       DO 33 JJ=I,NM,M
  205.       R(JJ)=R(JJ)+PIV*R(J)
  206.    33 J=J+M
  207.    34 KST=KST+MC
  208.       IF(ILR-MR)36,35,35
  209.    35 IC=IC-1
  210.    36 ID=K-MR
  211.       IF(ID)38,38,37
  212.    37 KST=KST-ID
  213.    38 CONTINUE
  214. C     END OF DECOMPOSITION LOOP
  215. C
  216. C
  217. C     BACK SUBSTITUTION
  218.       IF(MC-1)46,46,39
  219.    39 IC=2
  220.       KST=MA+ML-MC+2
  221.       II=M
  222.       DO 45 I=2,M
  223.       KST=KST-MC
  224.       II=II-1
  225.       J=II-MR
  226.       IF(J)41,41,40
  227.    40 KST=KST+J
  228.    41 DO 43 J=II,NM,M
  229.       TB=R(J)
  230.       MZ=KST+IC-2
  231.       ID=J
  232.       DO 42 JJ=KST,MZ
  233.       ID=ID+1
  234.    42 TB=TB-A(JJ)*R(ID)
  235.    43 R(J)=TB
  236.       IF(IC-MC)44,45,45
  237.    44 IC=IC+1
  238.    45 CONTINUE
  239.    46 RETURN
  240. C
  241. C
  242. C     ERROR RETURN
  243.    47 IER=-1
  244.       RETURN
  245.       END
  246.